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Abstract. Popular high-order schemes with compact stencils for Computational Fluid Dynamics 
(CFD) include Discontinuous Galerkin (DG), Spectral Difference (SD), and Spectral Volume (SV) 
methods. The recently proposed Flux Reconstruction (FR) approach or Correction Procedure using 
Reconstruction (CPR) is based on a differential formulation and provides a unifying framework for these 
high-order schemes. Here we present a brief review of recent developments for the FR/CPR schemes as 
well as some pacing items. 


1. Introduction 

In the field of Computational Fluid Dynamics (CFD), low -order methods are generally robust and 
reliable; as a result, they are routinely employed in practical calculations. For the same computing cost, 
high-order methods can provide considerably more accurate solutions, but they are more complicated and 
less robust. The need to improve and develop new high-order methods with favorable properties has 
attracted the interest of many researchers as evidenced by the recently held First (2012) and Second 
(2013) International Workshops on High-Order CFD Methods. 

The Discontinuous Galerkin (DG) method is currently among the most widely used high-order 
numerical methods for solving the compressible Navier-Stokes equations on unstructured meshes. It was 
introduced for the neutron transport equation by Reed and Hill (1973), analyzed by LaSaint and Raviart 
(1974) and developed and made popular for fluid dynamics equations by Cockburn, Shu, Bassi, Rebay, 
and others (see e.g., Cockburn, Karniadakis, and Shu 2000, Bassi and Rebay 1997a, b, Cockburn and Shu 
2005, 2009, Shu 2012, and the references therein). Efficient DG schemes storing data at nodal points 
known as nodal DG methods can be found in (Hesthaven and Warburton 2008). 

Alternative approaches to high-order accuracy employing the differential form (as opposed to DG 
which employs the integral form) have been proposed. Kopriva and Kolias (1996) pioneered this 
approach with the staggered- grid method on quadrilateral meshes. It was extended to triangular meshes 
by Liu, Vinokur, and Wang (2006) and named spectral difference (SD). Solutions for a wide range of 
problems by these methods can be found in (Wang et al., 2007 and Liang et al. 2009a and 2009b). 
Another class of schemes called spectral volume (SV) presented by Wang, Zhang, and Liu (2004) is 
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based on the idea of subdividing each cell into subcells or control volumes in a structured manner. A 
review of these as well as other types of high-order schemes can be found in (Wang 2007). 

Recently, an approach to high-order accuracy with the advantage of simplicity and economy called 
flux reconstruction (FR) was introduced in (Huynh 2007, 2009a). The approach amounts to evaluating the 
derivative of a discontinuous piecewise polynomial function by employing its straightforward derivative 
estimate together with a correction, which accounts for the jumps at the interfaces. The FR framework 
unifies several existing schemes : with appropriate choices of correction terms, it recovers DG, SD, as well 
as SV, and the FR versions are generally simpler and more economical than the original versions. In 
addition, the approach results in numerous new methods that are stable and super accurate, i.e., more 
accurate than expected (also known as super convergent). It was applied to ordinary differential equations 
in (Huynh 2009b), and the result is the Radau IIA collocation method. 

Whereas extension to a quadrilateral mesh is straightforward via tensor product, that for unstructured 
triangular meshes is nontrivial since tensor product no longer applies. Here, Wang and Gao (2009) 
showed that the derivative correction can be extended without the reconstruction concept. The resulting 
method was applied to solve the 2D Euler equations and named LCP (lifting collocation penalty). 
Extension to the 2D Navier-Stokes equations on meshes of mixed elements was carried out in (Gao and 
Wang 2009) and (Gao, Wang, and Huynh 2013). Extension to the 3D Euler and Navier-Stokes equations 
on mixed meshes was presented in (Haga, Gao, and Wang 2010, 2011, and Wang, Gao and Haga 2011) 
and to dynamic meshes by Yu, Wang, and Hu (2012). It was shown in (Huynh 2011) that the 
reconstruction concept applies to triangles as well. Due to the tight connection between FR and LCP, the 
involved authors combined the names and call them the CPR method (Correction Procedure via 
Reconstruction). PnPm-CPR schemes for the Navier-Stokes equations were studied in (Shiet al. 2012). A 
modification to assure that the resulting method is conservative was presented by Gao and Wang (2013). 
Adjoint-based error estimation and hp-adaptation were carried out in (Shi and Wang 2013), and 
comparisons for various types of schemes versus CPR were discussed in (Yu and Wang 2013). 

A mathematical foundation for the FR approach was recently provided by Jameson (2010), who 
proved that a particular SD scheme (recovered via FR) is energy- stable for ID linear advection. Vincent, 
Castonguay, and Jameson (2011a) subsequently extended this result, and proved that a one-parameter 
family of FR methods is energy-stable for linear advection. This family, referred to as Energy Stable Flux 
Reconstruction (ESFR) or Vincent-Castonguay-Jameson-Huynh (VCJH) schemes, includes a nodal DG 
method, the SD scheme previously identified by Jameson as being energy-stable, and the ‘g 2 ’ FR scheme. 
In further theoretical studies, VCJHschemes were extended to linear advection problem on 2D triangular 
grids by Castonguay, Vincent and Jameson (2012), to linear-advection-diffusion problem in ID by 
Castonguay et al. (2013), and to linear-advection-diffusion problem on 2D triangular grids by Williams et 
al. (2011) and Williams et al. (2013). Performance of VCJH schemes was also investigated using von 
Neumann analysis by Vincent, Castonguay and Jameson (2011b), and issues of nonlinear stability were 
discussed by Jameson, Vincent and Castonguay (2012). In the latter, it was shown that the location of the 
solution points is critical in terms of controlling aliasing driven instabilities. 

Next, we discuss some recent additional contributions to the development of FR/CPR. Suppressing 
oscillations near shocks via localized artificial diffusivity or LAD was carried out in (Miyaji2011) and 
(Haga et al. 2013). In the latter, FR methods for body-fitted Cartesian unstructured grids were also 
developed. Extensions of a particularly simple FR/CPR scheme called g 2 to the Navier-Stokes equations 
on moving and deforming domains were presented by (Liang, Miyaji, and Zhang 2013). A comparison of 
computational efficiencies of SD and CPR methods were carried out in (Liang, Cox, Plesniak 2013). 
Applications to turbulent internal flows for turbomachinaries and mitigation of aliasing errors were 
discussed in (Lu, Yuan, and Dawes 2012 and Lu, Liu, and Dawes 2013). An interface element approach 
dealing with non-conforming polynomials together with p-adaptation for viscous compressible flow 
simulations was elaborated in (Cagnone and Nadarajah 2012, Cagnone, Vermeire, and Nadarajah 2013). 
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In addition, an implicit large eddy simulation (ILES) solver was developed for CPR scheme using a third- 
order singly diagonal implicit Runge-Kutta scheme by Vermeire, Cagnone, and Nadarajah (2013). 

In this paper, we present a brief review of recent developments for the FR/CPR schemes as well as 
some key pacing items. Basic concepts of the DG and FR methods and their relations for a simple 
integration problem are presented in Section 2. Section 3 deals with conservation laws including the 
diffusion equation. Section 4 discusses 2D and 3D extensions of FR methods. Stability proofs are 
sketched in Section 5. Section 6 contains representative numerical results. Additional recent and ongoing 
research is described in Section 7. Some pacing items and areas for future research are discussed in 
Section 8. Finally, conclusions and discussion are presented in Section 9. 

2. A Simple Integration Problem 

We present the key ideas of both the integral and differential formulations and their relations using the 
following simple integration problem. On / = [—1,1], solve 

u'(x) = fix), u(— 1) = u L . (2.1) 


The exact solution is 


“exact GO = “l + J /(< Odf- 

The standard DG formulation is concise and results in a matrix equation discussed later, but it does not 
convey some key properties of the solution. Here, we will show that the DG solution u h of degree K — 1 
can be obtained as follows. First, project / onto the space of polynomials of degree K — 1 resulting in f h \ 
next, set U (x) = u L + / /),(/ )d<f; thus, U is of degree K. The solution u h is determined by the values of 
U at the K right Radau points (more on these points later). 

Examples. The proof for the above DG solution will be carried out after the following two examples, 
which convey the behavior of solutions. For the first example, find the linear DG solution of 

u’(x) = 2x, u(— 1) = 1. (2.2) 

By the above solution procedure, the linear projection of / is f h = f = 2x. Next, U (x) = 1 + 
2<fd<f = x 2 . The two right Radau points are —1/3 and 1. The linear DG solution is u h = ^x 
shown in Fig. 2.1(a). 

The second example is: find the cubic DG solution for 

u'(x) = 4x 3 , u(— 1) = 1. (2.3) 

Here, the cubic projection of / is f h = / = 4x 3 ; next, U (x) = u exact (x) = x 4 . The cubic DG solution is 
determined by the values of U at the four right Radau points shown in Fig. 2. 1(b). 

In both examples, U provides an approximation one degree higher than u h , and the values for both at 
the K right Radau points are the same. For this integration problem, since the ‘wind direction’ is unique, 
from left to right due to u L , the polynomial U is well defined. For the case of the Euler equations, 
however, the ‘wind direction’ is not unique and the definition of U becomes unclear. As will be shown, 
the FR technique of constructing U can still be applied to the fluxes. Also note that the gain in accuracy at 
the Radau points for the DG method was studied in (Adjerid et al. 2002) in the context of error estimates. 
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Figs. 2.1 (a) Linear DG solution for problem (2.2) and (b) cubic DG solution for problem (2.3). 


Preparations. Projections and orthogonal polynomials play an important role. Therefore, they are 
briefly reviewed. With any nonnegative integer m, let P m be the space of polynomials of degree m or 
less. For any two functions v and w on / = [—1, 1], denote their inner product by 


O ,w) 


r 


v(x)w(x) 


dx. 


The L 2 norm of v is 




A polynomial v is orthogonal to P m on / if, for each /, 0 < l < m, 



v(x) x l dx = 0. 


The criterion of orthogonality to P m provides m + 1 conditions. 

Legendre Polynomials. For each integer m > 0, the Legendre polynomial P m on / is defined as the 
unique polynomial of degree m that is orthogonal to P m _ t and P m (l) = 1. The Legendre polynomials 
are given by a recurrence formula (e.g., Hildebrand 1987): 

P 0 (x) = 1, P 1 (x)=x, 

and, for k > 2, 

2k — 1 k-1 

PkO) = — ^ — x P fe -i(x) — P k _ 2 (x). 

The first few Legendre polynomials are 

P 0 (x) = 1, Pi(x) = x, P 2 (x) = | (3x 2 - 1), P 3 (x) = i(5x 3 - 3x), P 4 (x) = i(35x 4 - 30x 2 + 3). 
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They are plotted in Fig. 2.2(a). 

Useful properties of the Legendre polynomials are listed below. If k > m, then P k is orthogonal to 
P m . Next, P k is an even function (involving only even powers of x) for even k, and an odd function for 
odd k. For all k, 


P k {~ l) = (-l) fe , P k { 1)=1. 


In addition, 


\\P k \\ 2 = (P k ,P k ) = 


2k +1' 


For later use, the derivative values at the end points are 


p;(-i) = ( ~ 1)t -‘ w+1) . />„■( i)=Mf±F. 


(2.4) 


The zeros of P k are the k Gauss points on [—1, 1]. 

Denote the projection operator onto P m by Pr m . Then for any function w, 


Pr, 


in 

t O) = ^ 

fc =0 


0 w,P k ) 
IIPJ 2 k - 


(2.5) 


Legendre Polynomials 

U 


Right Radau Polynomials 

U 



Radau polynomials. The right Radau polynomial of degree k (k> 1) is defined by 

(-D k 


Pr, k ~ 


<P k ~P k - 1)- 


( 2 . 6 ) 
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The zeros of R R k are the k right Radau points on [—1, 1]. The first few right Radau polynomials are 
plotted in Fig. 2.2(b). The above definition implies that R Rk is orthogonal to P k - 2 since both P k and 
P k _ 1 have this property. In addition, 

Rr,k(~ 1) = 1 and R Rk ( 1) = 0. (2.7) 

Note that R R k , which is of degree k, is defined by the above two conditions together with the k — 1 
conditions that it is orthogonal to P fe _ 2 . We can think of R R k as a polynomial approximating the step 
down function s defined by s(— 1) = 1 and s(x) = 0 for — 1 < x < 1. At x = ±1, 

k 2 (— k 

Rr.r'C- D = -y and R R , k \ 1) = . (2-8) 

The above expression for R r k \— 1) is related to the time step size restriction (or CFL condition) 
proportional to 1/k 2 . 

Lobatto polynomials. The Lobatto polynomial of degree k (k > 2), is defined by 

Lo k =P k -P k -2- (2-9) 

The zeros of the Lobatto polynomial of degree k are the k Lobatto points; they include the two boundaries 
+ 1 . The first few Lobatto polynomials are shown in Fig. 2.3(a). 


Lobatto Polynomials 



(a) 


Lagrange Polynomials 

u 



Figs. 2.3 (a) Lobatto polynomials and (b) Lagrange polynomials for the four Lobatto points. 


Lagrange polynomials. Let K be an integer > 1. For the nodal methods (e.g., Hesthaven and 
Warburton 2008), denote the nodes or solution points on / by k = 1 , K, which are typically the 
Gauss or Lobatto points. Let l k be the associated Lagrange polynomials, i.e., for each k, l k is of degree 
K — 1 taking on value 1 at ^ k and 0 at all other <f m , 


■ c °-n£fe 


m- 
m^k 


(2.10) 
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Lagrange basis functions for the case of four Lobatto points are shown in Fig. 2.3(b). 

DG formulation. Let the basis functions <p k for P K _ t be either the Lagrange polynomials l k (nodal 
DG) or the Legendre polynomials (modal DG) P k _ l5 k = 1, K. The solution is approximated by 


u h = 


K 

^ u k 4>k- 

k = 1 


( 2 . 11 ) 


Let <p be a test function, which is in P K -x ; typically, (f) is one of the basis functions <p k , k = 1, ... , K. 
Then the weak form for (2. 1) or it' = / is, formally (the correct equation is (2. 16) or (2. 17)), 

(u',0) = (/,0). (2.12) 

Let f h be the projection of / onto P K _ f h = Pi V-i(/) via (2-5) with m = K — 1. Since <p varies on 
P K ~ i, wehave ( f,<p ) = (f fl , 0), i.e., we can project first, then apply the inner product. The above implies 

(u',0) = (/ h ,0). (2.13) 

Going back to (2. 1), the solution of u' = f h and it (— 1) = u L is 

U(x) = u L + J f h (Odf- (2.14) 


Clearly, U is of degree K since f h is of degree K — 1. But the solution we are seeking is u h of degree 
K — 1. The question is: how do we obtain u h from U1 Or, to put it differently, how does the DG solution 
u h relate to U1 

To answer the above question, we first define the DG solution. For the DG method, u h {— 1) is allowed 
to be and is generally different from u L . Applying integration by parts to the left hand side of (2.13), 

[“/i0]-i - (u h ,0') = (/ h ,0). (2.15) 

To involve the initial condition, we replace it h (— 1) in the boundary term above by u L . Thus, the DG 
solution u h is required to satisfy, for all 0 in P K _ l5 

u h (l)0(l)-u L 0(-l)-(u h ,0') = (/ h ,0). (2.16) 

Applying integration by parts again to (u h ,<p'), 

~[u L - u h (- 1)] 0(- 1) + ( u h ',(t) ) = (/ h ,$). (2. 17) 

This time, for the boundary term, it h (— 1) is employed. The above, the result of integrating by parts twice, 
is the ‘strong form’ whereas (2.16) is the ‘weak form’ (Hesthaven and Warburton 2008). The solution can 
be obtained by solving either (2.16) or (2.17) where u h is of the form (2.11), and <p is replaced by <p k , 
k = 1, ...,K\ the result is a matrix equation for the K unknowns u k . 

Flux Reconstruction formulation. Our goal is to eliminate the test function in (2.17) so that the 
integral formulation results in a differential one. To this end, we raise the following question: can we find 
a polynomial g LB = g on / which possesses the property that for any (f) of degree K — 1 or less, 


-</>(-!) = (g',4>)- 


(2.18) 
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If such a polynomial exists, we can combine (u' h , 0) in (2. 17) with (g',0) and eliminate the test function. 
From (2.17), since f h is of degree K — 1, we require g' to also be of degree K — 1; as a result, g is of 
degree K. Applying integration by parts to the right hand side of (2.18), we have 


-</>(-!) = g( 1)0(1) - s r (-i) 0 (-i) - ( g , 0 ')- 


The above holds if g satisfies, 


#(-!)= 1, g( 1) = 0, 


(2.19) 


( 2 . 20 ) 


and, for all 0 in P K _ l5 


(<7,0')=O. (2.21) 

Since 0 is of degree K — 1, 0' is of degree ft' — 2; moreover, 0' spans P K - 2 as 0 spans P K _ 1 . The above 
condition implies that g is orthogonal to P K - 2 , i.e., for any polynomial cp of degree K — 2, 

(g,(p) = 0 . ( 2 . 22 ) 

Orthogonality to P K _ 2 provides K — 1 conditions; (2.20) provides the other two. Loosely put, the 
condition g(— 1) = 1 deals with the jump at the left boundary whereas the above together with g(l) = 0 
implies that away from the left boundary, g approximates the zero function. By the discussion after (2.7), 
g is the right Radau polynomial, 


9 = Rr,k ■ (2.23) 

Thus, the answer to the question posed for (2. 18) is positive. Note that g = g LB , which corrects for the 
jump at the left boundary, is the right Radau polynomial (vanishing at x = 1). 

We now return to our task of eliminating 0. With g as above, ( g' ,<fi ) = — 0 ( — l) - Thus, by (2.17) 

[u L -u h (- 1)] (g',<p) + (u00) = (/ h ,0). (2.24) 

What is crucial here is that 0 can be factored out and canceled: 

([u L -u h (-l)] 5 r + u h )'= f h . (2.25) 


Therefore, [u L — u h (—l)]g + u h is a polynomial of degree K whose value at the left boundary is 

[u L - u h (-l)]g(-l) + u h {- 1) = u L . (2.26) 

The above two conditions imply, with U defined by (2.14), 


[u L -u h (-l)]g + u h = U. 


(2.27) 


Next, denoting the K right Radau points by r k , k = 1, ..., K. Since g vanishes at these points, 


u h (r k ) = U(r k ). (2.28) 

Thus, with U of degree K defined by (2. 14), the DG solution u h of degree K — 1 is defined by the values 
U at the K right Radau points. 
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The key idea of FR is the following. Given u h which does not match the boundary condition u L , we 
can add the correction term [u L — u h (— 1)]$ to u h so that the resulting U takes on the left boundary value 
u L , does not alter the value u tl (l), and retains the property of ‘best possible’ approximation to u h . The 
function U is called the ‘continuous flux function’ as opposed to u h , which is discontinuous across the 
cell boundaries. It turns out that for conservation laws, there are numerous ways to define g so that the 
resulting scheme is stable. 

For the ordinary differential equation (ODE) it'(x) = f(x, u), u(0) = u 0 , there is only one ‘wind 
direction’, i.e., x moves forward. The above argument shows that for such a case, it is sensible to use the 
K right Radau points as collocation or quadrature points, which implies making use of U. The resulting 
scheme is identical to an implicit collocation Runge-Kutta scheme called Radau IIA (Huynh 2009b). 

3. FR/CPR Methods for the One-Dimensional Case 
Conservation laws. Consider the conservation law 


u t +f x = 0 (3.1) 

with initial condition u(x, 0) = u init (x) and the flux / depends on it. The solution u is assumed to be 
periodic or of compact support so that boundary conditions are trivial. 

Let the domain of calculation be divided into (possibly nonuniform) cells or elements E p j = 1, 2 , ... 
Denote the center of Ej by Xj and its width by hj. With <f varying on / = [— 1, 1] and x on Ej, the linear 
function mapping / onto Ej and its inverse are 

x(<f) = Xj + %hj/2 and ^ (x) = 2 (x — Xj)/hj . 

In addition, denote the nodes or solution points on /, which are typically the Gauss or Lobatto points, by 
<f fe , k = 1, K. They relate to the nodes on Ej by 

x j,k ~ x (£k) ~ x j "b ^k^j /2. (3.2) 

The global derivative, e.g., f x , can be obtained from the local one via the chain rule — = ~~~ ■ 

dx hj dq 

Whereas the FR approach can be formulated in modal form, for simplicity, we discuss only the nodal 
form. On Ej, with the Lagrange basis functions <fij,k( x (.0) = k = 1, ..., K, let the solution be 

approximated by a polynomial of degree K — 1 denoted by u 

K 

u j = ^ u j,k4>j,k- ( 3 - 3 ) 

fc = i 


Here, Uj k depends on t, and cf>j k , on x. 

At time level n, suppose Uj k are known for all j and k. We wish to calculate dUj k /dt = 
dUj k (t) / dt at t = t n , i.e., to calculate {f x ) j, k - Then, we march in time by, say, a Runge-Kutta method. 

With fj k = f(Uj k ), let fj (x) be the polynomial of degree K — 1 interpolating fj k , k = 1, ... , K, 
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(3.4) 


K 

fj = X fj, k 4>j,k- 

/e= i 

The flux polynomials {fj} form a function, which is generally discontinuous across cell interfaces and fj 
is called the discontinuous flux function. Note that ( fj ) involves no interaction of data among cells. 

To account for interaction, we construct a continuous flux function, which approximates the 
discontinuous function in some sense, and then calculate its derivative. The continuous flux function will 
be obtained by adding a correction to the discontinuous one. As a result, we still need the derivative of the 
discontinuous function. At each <f fc , k = 1 , ... ,K, it is easy to derive the derivative matrix {d kl } where 

K 

= 0 - 5 ) 

fc= 1 

Instead of differentiating fj as above, we can use the chain rule 

(/f) yic = a i u j. k )i u ^ jk - (3.6) 

Wang and Gao (2009) found that for the Euler equations, the chain rule yields more accurate solutions. 

At each interface x y+1 / 2 , set 

u l = u j+i /2 ,L = tty(l) and u R = Uj +1/2 ,R = Uj +1 (-l). (3.7) 

From these two values, we can obtain a common flux (shared by the two adjacent cells) denoted by 
f) +i/ 2 , com- For advection problems, the common flux is typically an upwind flux; for diffusion problems, 
however, it is usually a centered quantity. This flux is often called the ‘numerical flux’ by the DG 
community. 

Next, we reconstmct the flux by a continuous function F such that on each cell Ej, F is a polynomial 
denoted by Fj approximating the discontinuous flux function fj. To assure continuity across cells, Fj is 
required to take on the common flux values at the two interfaces: 

Fj \Xj— 1 / 2 ) fj-i/2, com and Fj(jXj+i/2j fj+ 1 / 2 , com- (3.8) 

In addition, Fj is required to be of degree K so that its derivative is of degree K — 1, matching that of u h . 
Switching to the local description, 

Pji-f)~ fji~f) = fj-i/2, com -//(-!) and Jy(l) - /y(l) = f J+1/2i com - / ; -(l). (3.9) 

Therefore, Fj — fj takes on the above prescribed left and right correction values, is of degree K, and 
approximates the zero function in some sense. 

We now separate the prescription of the correction at the left interface from that of the right. On /, let 
g LB be the correction function for the left boundary defined by 

g LB (- 1)=1, fl LB (D=0, (3.10) 
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and g LB is a polynomial of degree K approximating the zero function in some sense. Let g RB be the 
correction function for the right boundary defined by reflection 

9rb(0 = 9lb(~' ()■ (3.11) 

For the left interface Xj_ 1 / 2 , the polynomial 

fj(0+ [fj -1/2, com fj ( 1)] #lb(0- (3.12) 

provides a correction to fj (<f) by changing the flux value at this interface from fj (—1) to fj- 1 / 2 , com while 
leaving the value at the right interface unchanged, namely fjiX). Next, the polynomial 

Fj ) = /;■ (O + [fj- 1/2, com - fj (-1)] 9lb (0 + [fj+ 1/2, com “ fj (1)] 9rb (0 (3. 13) 

provides corrections to both interfaces: F ; (— 1) = f_ i/ 2 , C om an d F ; (l) = //+ 1 / 2 , com- Thus, Fj is of 
degree K, takes on the two common flux values, and approximates fj in the same sense that g LB and g RB 
approximate the zero function. The derivative of Fj at the solution point / k is 

( F /) jk = ( fdj, k + com -/;(-!)] WCffc) + [/;■+ 1 / 2 , com - /yCO] 9rb' ( fk) (3.14) 

The derivative ( F x )j k follows. The solution Uj k can then be updated via, say, a Runge-Kutta method. 
What is crucial in (3.14) is that at each solution point, the derivative (/L) . of the continuous flux 
function is obtained by correcting the derivative (/^) . of the discontinuous flux function. The correction 

amount is straightforward once the values g LB and g RB ' (f k ) are known. These derivative values, in 
turn, can easily be derived once g LB and g RB are defined on I. 

We summarize the FR/CPR algorithm below. 

Algorithm. At time level n, suppose Uj k are known for all j and k. 

(1) At each interface j + 1/2, if the left and right values of u are not available, calculate them; then 
estimate and store the common (upwind) fluxes at all interfaces. 

(2) In the cell j, for k = 1, K, evaluate fj k = f(uj k ); then obtain (/?) . of the discontinuous 
flux function by (3.5). Alternatively, the chain rule (3.6) can be employed. 

(3) At the two interfaces of Ej, get the corrections fj- 1 / 2 , com — fj ( — l) cmd fj+ 1 / 2 , com — fj (!)■ At the 
solution points, evaluate (Ff). by (3.14) and then (F x )j k . 

(4) March in time by, say, a Runge-Kutta method. This completes the algorithm. 

Correction functions. Next, we discuss various choices for the correction function g = g LB . Recall 
that g is determined by (3. 10) together with the K — 1 conditions. These different choices for g result in 
the DG, SD, SV as well as new schemes. Thus, the FR approach provides a unifying framework for these 
methods. 

The first choice for g requires that it is orthogonal to P K - 2 , which means it is the right Radau 
polynomial as in (2.23), and the resulting scheme is identical to DG, 
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9 — 9t>g — Rr.k- 


( 3 . 15 ) 


The condition of orthogonality to P K - 2 can be relaxed. It was verified via Fourier analysis (Huynh 
2007, 2009a) that if g is orthogonal to P K - 3 the resulting scheme is stable (the converse is not true, 
however). Since both R r k and R r K _ l are orthogonal to P K - 2 , such a correction function can be written 
as, 


9 — aR R K + (1 cc)R r k _ 1 (3.16) 

where, 0 < a < 1 and a remains to be determined. 

The second choice for g, denoted by g 2 or STump, lo (for Tumping for Lobatto points ’), is defined as 
follows. Since a steeper correction function tends to result in a scheme with a smaller CFL limit, to make 
g less steep, the extra condition is obtained by pushing one of the zeros to the right boundary, Le., f = 1 
is a zero of multiplicity two. After some algebra, 


9 2 STurnp, Lo 


K-l K 

2K^1 Rr ’ k + 2K^1 Rr ’ k ~ 1 - 


(3.17) 


The function g 2 has the following remarkable property. Among the K Lobatto points, g 2 vanishes at 
K — 1 of them; the exception is the left boundary. 

The final choice here for g requires that g van is hes at the K — 1 Gauss points: 


9 Ga — 


K K-l 


(3.18) 


Whereas the staggered- grid and SD schemes are mildly unstable, the above provides a modification using 
the K — 1 Gauss points as flux points. The resulting scheme is stable for all K. 

It can be shown by calculations using Fourier analysis that the scheme using g DG is stable and accurate 
to order 2 K — 1, and the schemes using g 2 and g Ga are stable and accurate to order 2 K — 2. In general, if 
g is orthogonal to P m , the resulting scheme is accurate to order K + 1 + m. If g is orthogonal to P K - 2 , 
i.e., g is of the form (3. 16), the resulting scheme is Fourier stable. Stability proofs for FR schemes will be 
discussed in Section 5. 

Also note that the steepest slope of g, which often takes place at the left boundary, relates to the time 
step size limit. For example, for DG, as in (2.8), g' = R' r fe (— 1) = —k 2 / 2; when an explicit Runge-Kutta 
method is employed, it is well known that the time step size limit for the DG scheme of degree /c — 1 is 
roughly proportional to 1/k 2 . 

The Diffusion Equation. On (— oo, oo), consider the diffusion equation, 

Uf u xx (3.19) 


with initial condition 


it(x, 0) = u 0 (x) (3.20) 

As in the case of (3.1), assume that the data Uj k are known at time level n. We wish to evaluate the 
second derivative in a manner which takes into account the data interaction among cells. For simplicity 
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and efficiency, the stencil of the scheme is required to remain compact in the sense that the second 
derivative evaluation in a cell involves the data of only that cell and the two immediate neighbors. 

Common values and corrected derivative estimates. The first task is to estimate u x at the solution 
points Xj k . Since the function {it ; } is discontinuous across the interfaces, to estimate u x , we first 
reconstruct it by a piecewise polynomial function {uj }, which is continuous across the cell interfaces , and 
on each Ej, is of degree K and approximates it ; - (the superscript ‘C’ stands for ‘continuous’ or 
‘corrected’). The derivative approximation (it| : ) ( Xy fe ) accounts for the data interaction. 

In order for {itj'} to be continuous at the interfaces, u c j and uj +1 must take on the same value at x J+1 / 2 . 
Thus, at each interface, we need to define a common interface value (or common value). Here, for a 
diffusion problem, we use a centered-type quantity: with u L and u R given by (3.7), 

u com~ u y+l/ 2 , com — ( U L T U^)/2. (3.21) 

The above formula was employed by Bassi and Rebay (1997a, 2000). A more general formula is the 
weighted average, with 0 < k < 1, 

u com — u j+ 1 / 2 , com — K U L T (1 — k)u r . (3.22) 

For k = 0 or k = 1, we have the one-sided formula used in the local DG or LDG (Cockburn and Shu 
1998) as well as the compact DG or CDG methods (Peraire and Persson 2008). 

Next, we require uj (x) to take on the common values Uj _ 1 / 2 com at x j- 1/2 ar| d Uj+ 1 / 2 , com at x j+ 1 / 2 > 
to be of degree K, and to approximate n ; (x). That is, in the local coordinate, 

Uj (0 — Uj (0 + [if/ — 1 / 2 , com 1)] 0LB (?) T [ 19 + 1 / 2 , com dRB (O- (3.23) 

The derivative is 

= (u/) f (0 + [u y _i/2,com- “;(-!)] 9lb' (. 0 + Wj+ 1 / 2 , com - UyCD] 9rb' (O- (3.24) 
And the derivative (itf) follows. See Fig. 3.1(a). 

J X 



0 12 3 



0 12 3 


( a ) 


(b) 


Figs. 3.1. Centered-type common derivative: (a) using a four-cell stencil and (b) using a two-cell 
stencil via (3.27). Here, the solution polynomials are linear, and the correction function g DG is parabolic. 
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Common derivative and corrected second derivative estimates. At each interface, in formula 
(3.21) for the common value, with 0 < k < 1, the weight for u L is k and that for u R is 1 — k. To define 
the common derivative value, we switch the two weights. Loosely put, this switch makes the method 
unbiased. If we apply the weighted average to (uj) , the resulting (u x )j+i/ 2 , com has a stencil of four 
cells, from j — 1 to j + 2 (see Fig. 3.1(a)). Since the calculation of u xx in cell j employs (u x )j_ 1 / 2 , CO m 
and (Ux)j+i/ 2 , com’ the corresponding scheme has a five-cell stencil. 

We now define a common derivative at j + 1/2 that involves only the data in the two adjacent cells. A 
scheme with such a compact stencil is desirable since it is easy to code, the boundary conditions involved are 
simple, and the resulting implicit version has a sparse and generally invertible matrix. To this end, correcting 
for the right boundary of cell j, set 

uf B (0 = Uj(0 + [Uj +1 /2,com-Uj(l)]g RB (0 (3.25) 

i.e. , uf B corrects for the right boundary, namely uJ B (1) = Uj+ 1 / 2 , com- while leaving the value at the left 
boundary unchanged, namely, Uj(— 1). Next, correcting for the left boundary of cell j + 1, set 

u j+i (C) — u j+i(0+ [wy+i^, com u j+i( 1)]/7lb( < D (3.26) 


Then corrects for the left boundary, uj^ (— 1) = tty +1 / 2 ,com- while leaving the value at the right 
boundary unchanged, namely Uy(l). 

Finally, for the common derivative at j + 1 /2, set 


( u x ) 


xJj+1/2, com 


+ 


2 c 

u j+ 1/2, com 11 j (1)] S'rb'CI)) 

K [“;'+i/2, com -“;■+! (-1)]^lb'(-1)}- 


(3.27) 


See Fig. 3.1(b). Note the dependence only on tt ;+1 / 2 ,com an d the data on Ej and Ej+ 1- 

With the corrected derivative given by (3.24) and the common derivative above, we can obtain the 
corrected second derivative estimates. 

The above procedure yields the CPR versions of the BR2 scheme if k = 1/2 (Bassi and Rebay, 
2000) and the LDG (Cockburn and Shu 1998) or CDG schemes (Peraire and Persson 2008) if k = 0 or 
K = 1. Comparing the centered versus one-sided common values with g DG as correction function (i.e., 
BR2 versus CDG), the former is of order 2 K — 1; the latter, order 2K\ the former, however, has the 
advantage that its CFL limit is more than two times larger than the latter. Since super -convergence (or 
super- accuracy) does not hold for the general case of nonlinear equations, the scheme using centered 
common values appear to have a slight edge. 

4 Two/ Three-Dimensional Extension 


The extension of the CPR formulation to quadrilateral and hexahedral elements is straightforward. The 
basic idea is to first transform the governing equations from a physical element to the reference or 
standard element. Then, the ID CPR formulation is applied on the standard element in each coordinate 
direction. 

Consider the 2D conservation law 
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u t + fx + 9 y — 0 - 


( 4 . 1 ) 


Denote f = (x, y) the coordinates of the physical domain, and <f = (£, 77 ) the coordinates of the standard 
element. Let the transformation be written as 


r = 


X «, (f) 




(4.2) 


Where ry are the physical coordinates used to define an element, and Mj (<f ) is the shape function. The 
transformed equation takes the following form 


u t + + g v — 0 (4.3) 

where ft = |/|.it, /= I/I ■ {% x f + <7 = I/I ■ (77*/ + Vyd) and / is the Jacobian matrix of the 

transformation, i.e., 



x v 

y* ■ 


(4.3) 


For a quadrilateral element of index / (not related to / of (4.2)), two indices (k,m) are used to denote 
the solution point, and Hi. km denotes the degrees of freedom (DOFs). The CPR formulation is then 

duj. k m dfj. k m dg j. k rn 
dt + + dg 

+ [/com (-!<%, m) -fj (-!<%, m)]^ + \fcom{^9k,m) ~ fj (Mte.m)] (4 ‘ 4) 

+[gcom(Sk,m'- 1 )-9j(tk,m'- 1 )]-^ + [/com( $k,m > - 1 ) “ 9j( Sk.rrv - 1 )] 

where the constants <x L and a R are the derivatives of the correction functions, and are also called 
correction coefficients. 

The extension to simplex and other types of elements is not as straightforward since the correction 
functions are not readily available. The first extension to triangular elements was based on the so-called 
lifting collocation penalty approach (LCP) (Wang & Gao, 2009). As it turns out, the final form is very 
similar to (4.4). The details of the derivation are omitted and we summarize the basic formulation here. 
Define two sets of points, solution points and flux points as shown in Figure 4.1. The CPR formulation 
can be rewritten as 


du jk ^ dfjuj) ^ dg(Uj) 
8t + 3 x m + dy jjc 


1 



Z2> 

fedVj l 


k,f 


,[/"], A =°- 


(4.5) 


where |V ; | is the area of the triangle, and S f is the length of side/, and [/"] is the normal flux 

difference between the common Riemann flux and the internal flux. The extension to 3D elements 
follows a similar path, and can be found in (Haga, Gao, Wang 2011). 
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Figure 4.1. Solution points (squares) and flux points (circles) for a degree 2 element 

5 Stability of Flux Reconstruction Methods 
5.1 Energy Stability for Linear Advection Problems in ID 

Jameson (2010) recently proved that a particular SD scheme (recovered via FR) is energy-stable for 
linear advection problems in ID. Vincent, Castonguay, and Jameson (2011a) subsequently extended this 
result, and identified a family of stable FR schemes for linear advection problems in ID (for all orders of 
accuracy). Specifically, it was proven that if the left and right flux correction functions are defined as 
follows 


9l = 




L k + 


( r lk^‘k-l + ^k+l \ 

V 1 +Vk /J' 

fVk^k-l + ^fc+l\ 
l l+rj k )\' 


(5.1) 

(5.2) 


where k = K — 1 is the degree of the solution polynomial within each element, L k _ t , L k , and L k+1 are 
Legendre polynomials of the denoted degree (normalized such that |L fe (±l)| = 1 for all k), and 

c(2k + 1 )(a fe /c!) 2 
9k= 2 ' 

with 

(2/c) ! 

Uk ~ 2 fe (/c!) 2 ' 

and c a free parameter in the range 

(2/c + l)(a k k\) 2 <c<co ’ 

then a broken Sobolev type norm of the approximate solution is guaranteed to be non-increasing, and thus 
bounded. Consequently, by equivalence of norms in the finite-dimensional solution space, any norm of 
the solution is guaranteed to remain bounded, and thus the method is guaranteed to be stable. 

The resulting one parameter family of FR schemes, defined in terms of the free parameter c, have been 
referred to as Energy Stable Flux Reconstruction (ESFR) schemes or Vincent-Castonguay-Jameson- 
Huynh (VCJH) schemes. It can be noted that judicious choice of the parameter c leads to recovery of 
various known FR schemes. Specifically, if 

c = 0, 

then a nodal DG scheme is recovered, if 
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2k 

(2/c+ l)(/c + l)(a fe /c!) 2 ’ 

then a particular SD scheme is recovered (the scheme is, in fact, the particular SD scheme that Huynh 
showed to be Fourier stable, and Jameson (2010) proved to be energy stable), and if 

2{k + 1) 

(2k+ l)k(a k k\) 2 ’ 

then the g 2 FR scheme is recovered. 

5.2 Energy Stability for Linear Advection Problems on 2D Triangular Grids 

VCJH schemes for linear advection problems in ID have been extended by Castonguay, Vincent, and 
Jameson (2012) to treat linear advection problems on 2D triangular grids. As in the ID case, a one- 
parameter family of correction functions were identified that guarantee a particular norm of the solution is 
non-increasing. However, unlike in the ID case, an explicit expression for these correction functions was 
not presented (instead the divergence of each correction function was defined implicitly via a matrix 
system). Interestingly, the one-parameter family of schemes did not appear to include a SD scheme as a 
special case — despite the fact that Balan, May, and Schoberl (2012) were able to identify stable SD 
schemes on triangular grids for several orders or accuracy. 

5.3 Energy Stability for Linear Advection -Diffusion Problems in ID and on 2D Triangular Grids 

Recently Williams et al. (2011), Castonguay et al. (2013) and Williams et al. (2013) have extended 
VCJH schemes for linear advection problems to develop a range of VCJH schemes for linear advection- 
diffusion problems. Their approach involves use of VCJH correction functions to construct a C° 
continuous polynomial representation of the solution (in addition to a C° continuous representation of the 
flux) within each element. Development of an energy-stable treatment for diffusive terms is important, 
since it is a prerequisite for effective solution of the Navier-Stokes equations. 

5.4 Non-Linear Stability 

Jameson, Vincent and Castonguay (2012) showed that FR methods can be afflicted by an aliasing 
driven instability if the flux function is non-linear. Such instabilities are a consequence of aliasing errors 
(that occur when a polynomial representation of the non-linear flux is constructed via a collocation 
projection at the solution points). Jameson, Vincent and Castonguay (2012) also demonstrated that the 
location of the solution points plays a critical role in determining the extent of any aliasing driven 
instabilities. Specifically, they suggest that the solution points should be located at the abscissa of a strong 
quadrature rule in order to minimize aliasing driven instabilities. This finding is supported by the 
numerical experiments of Castonguay, Vincent and Jameson (201 1), who used the FR approach to solve 
the Euler equations on 2D triangular grids. They found that if the solution points were located at the so- 
called alpha- optimized points of Hesthaven and Warburton (2008) then the simulations blew up. 
However, if the solution points were located at the abscissa of a high-strength quadrature rule derived by 
Taylor, Wingate and Bos (2005), then the simulations remained stable (see Fig. 5.1). 

5.5) Von Neumann Analysis 

Energy-based stability proofs are powerful since they apply for all orders of accuracy and on non- 
uniform grids. However, they do not offer insight into all the stability properties of a numerical scheme. 
Huynh (2007) and Vincent, Castonguay and Jameson (2011b) presented comprehensive von Neumann 
analyses of FR methods in order to elucidate further stability properties of the schemes. Their results 
indicate that the form of the flux correction function has a significant impact on the CFL stability limit 
associated with a given FR scheme. In the context of ID VCJH schemes for linear advection, it has been 
shown that increasing the free parameter c (from zero) can increase the CFL limit by over a factor of two 
in certain cases (at the cost of a reduction in the overall accuracy of the scheme). 
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-4 -2 0 -4-2 0 


(a) (b) 

Fig.5.1 Density color maps from inviscid simulation of Euler vortex motion in a cross -flow. In both (a) 
and (b) the domain was meshed with triangular elements, and fourth-order solution polynomials were 
used to represent the solution within each element. However, in (a) solution points are located at the 
alpha-optimized points of Hesthaven and Warburton (2008) and in (b) solution points were located at the 
abscissa of a high-strength quadrature rule derived by Taylor, Wingate and Bos (2005). The solution in 
(a) quickly becomes unstable and the simulation blows -up, whereas the solution in (b) remains stable. 
Adapted from study of Castonguay, Vincent and Jameson (201 1). Copyright P. Castonguay, P. E. Vincent 
and A. Jameson. Reproduced with permission. 


5.6 Remaining Stability Problems 

The above developments are significant in terms of understanding fundamental stability properties of 
FR schemes. However, there remain various stability issues that need addressing. Firstly, whilst empirical 
evidence suggests that the ID stability proofs of Vincent, Castonguay and Jameson (2011a) and 
Castonguay et al. (2013) extend to tensor product elements (quadrilaterals and hexahedra), there exists no 
mathematical proof, and it remains an open question as to whether multi-dimension tensor product 
formulations based on ID VCJH schemes are in fact linearly stable. Moreover, as yet VCJH schemes for 
advection and advection-diffusion problems have not been extended to tetrahedral, prismatic or pyramid 
shaped elements, all of which are widely used to create unstructured meshes of complex 3D geometries. 
Finally, robust strategies for reducing/controlling aliasing driven instabilities in multiple element types 
need to be developed. 


6 Representative Numerical Examples 

In this section, two numerical examples are shown to demonstrate the capability of the CPR 
formulation. 

6.1 Direct Numerical Simulation of the Taylor-Green Vortex at Re = 1600 

This is a benchmark case (C3.5) of the 1st International Workshop on High-Order CFD Methods 
(Wang et al. 2013), and was designed to evaluate numerical methods in accurately capturing the evolution 
of a smooth flow to a turbulent flow. Avery high -resolution simulation with a spectral method is used as 
the “analytical solution”. In the present simulations, p2, p3 and p4 CPR schemes were used with different 
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mesh sizes ranging from 64 3 to 96 . Table 6.1 summarizes all the cases with the time steps used in the 
simulations. The 3 ld -order SSP (strong-stability preserving) Runge-Kutta scheme was used for time 
integration. The computed energy dissipation rate and enstrophy history are displayed in Figure 6.1. The 
obvious trend is that the higher-order schemes perform better than the lower order schemes. For example, 
the p3 scheme with 64 3 cells performs better than the p2 scheme with 96 3 elements, and is less expensive. 


Table 6.1 Summary of Cases 



Grid (Hex) 

P 

nDOFs 

Time step 

comp. 1 

64x64x64 

2 

7,077,888 

3.92e-04 

comp. 2 

64x64x64 

3 

16,777,216 

3.92e-04 

comp. 3 

64x64x64 

4 

32,768,000 

2.63e-04 

comp. 4 

96x96x96 

2 

23,887,872 

3.92e-04 



Figure 6.1. Computed histories of energy dissipation rate and enstrophy 
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6.2 Computations of bio-inspired vortex-dominated flows 

This case was performed by Yu et al (2012), and re-produced here with permission. Flows over a 
rectangular flapping wing, as shown in Fig. 6.2, are studied here. The wing can undergo a flapping or 
combined flapping-pitching motion. A remeshing technique is then used to deform the mesh at each time 
step. 



Figure 6.2. Wing surface and root plane meshes for a rectangular wing. 



Figure 6.3. Comparison of the vortex topology for the rectangular and bio-inspired wings at four phases 
(0°, 90°, 180° and 270°) with the flapping motion. The upper row is for the rectangular wing whereas 
the lower row is for the bio-inspired wing. 


In this study, the Strouhal number (St) of the finite-span flapping wing was selected to be well within 
the optimal range usually used by flying insects, birds, and fish (i.e., 0.2 < St < 0.4). The Mach number of 
the free stream is set to be 0.05 to mimic incompressible flow. The Reynolds number (Re) based on the 
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free stream velocity and the maximum chord length is 1200. The reduced frequency of the flapping 
motion is 3.5, and the Strouhal number of the wingtip is 0.38. The space discretization accuracy for the 
simulation is of third order, and the time integration is performed with the explicit third order TVD 
Runge-Kutta method. 

The flapping alone and the combined flapping-pitching were studied. The computed iso-surfaces of the 
Q-criterion at different phase angles for both motions are shown in Figure 6.3. The histories of the thrust 
coefficient from both motions are presented in Figure 6.4. Note that the thrust from the combined motion 
is more than an order of magnitude higher than that from flapping alone. In fact, the averaged thrust from 
the combined motion is 27 times as large as that from the flapping motion. 



Fig. 6.4 Thrust coefficient histories for the wing with the combined and the root-fixed plunging motions. 


7 Additional Recent Research on FR/CPR Methods 

In this section, we briefly discuss additional recent contributions to the development of FR/CPR. 

Suppressing oscillations near shocks via localized artificial diffusivity or LAD was carried out in 
(Miyaji 2011) and (Haga et al. 2013). Initial results are encouraging; however, dealing with shocks 
requires further research and tests, especially for the case of multiple dimensions. 

Concerning SD-type schemes, extensions of a particularly simple FR/CPR scheme called g 2 to the 
Navier-Stokes equations on moving and deforming domains were presented by (Liang, Cox, Plesniak 
2013) and (Liang, Miyaji, and Zhang 2013). They demonstrated that the g 2 scheme is (up to 40%) 
faster and easier to implement than the SD method. Comparison of computational efficiencies for 
various types of schemes versus CPR was carried out (Yu and Wang 2013). Their results show that CPR 
schemes are the most efficient. 

An interface element approach dealing with non -conforming polynomials together with p-adaptation 
for viscous compressible flow simulations was elaborated in (Cagnone and Nadarajah 2012, Cagnone, 
Vermeire, and Nadarajah 2013). In addition, an implicit large eddy simulation (ILES) solver was 
developed for CPR schemes using a third-order singly diagonal implicit Runge-Kutta scheme by 
Vermeire, Cagnone, and Nadarajah (2013). Preliminary work on FR methods for body-fitted Cartesian 
unstructured grids was presented in (Haga et al. 2013). A PnPm-CPR method for the Navier-Stokes 
equations was studied in (Shi et al. 2012). Here, among the main findings in applying PnPm to CPR was 
that the accuracy gain in multiple dimensions with a compact reconstruction stencil is rather limited. 
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It is well-known that hp-adaptation is beneficial for high-order methods (Wang et al. 2013). Adjoint- 
based error estimation and hp-adaptation for CPR was studied in (Shi and Wang 2013). Their results 
show significant savings compared to the uniform h or p refinement. 

Progress toward applying FR/CPR methods to solve practical internal flow problems has been made 
by Lu, Dawes, and Yuan (2012), Lu, Liu, and Dawes (2013), and Lu and Dawes (2013). These authors 
have devised a parallel solver on hybrid unstructured meshes including tetrahedra, prisms, pyramids and 
hexahedra for turbulent subsonic /transonic flows. 

8 Pacing Items 

Pacing items are similar to those mentioned by the committee for the International Workshop on High- 
Order CFD Methods (Wang et al. 2013). Some of the key items are discussed below. 

High-order mesh generation. In order to achieve high-order accuracy, curved geometries need to be 
represented with high-order polynomials. The generation of unstructured, highly-clustered viscous 
meshes near high-order boundaries requires further research to improve robustness. The main difficulty is 
that cells near the curved geometries can overlap each other. 

Capturing shocks. The two main approaches are local artificial dissipation (LAD) and limiting. The 
former involves user specified parameters, and latter often causes convergence to stall. A third approach 
is via h-p mesh adaptation where, near shocks, the mesh is refined and the method switches to the first- 
order upwind scheme; here, deciding when to switch is not trivial. An optimal method should capture 
shocks with high resolution, preserve accuracy at smooth regions of the flow, and be convergent when 
needed. 

Time stepping. Much research is needed on how to handle the stiffness generated by highly- 
anisotropic meshes near walls for viscous flows. Low storage and efficient iterative solution methods for 
both steady and unsteady flow problems are active area of research. Time stepping methods are derived 
typically by mathematics where stability and accuracy are the main focus. Can time-stepping methods be 
derived in combination with physics to deal with the stiffness issue? 

Mesh adaptation. There are numerous works on this topic, however, a simple, efficient, and robust 
3D adaptation method requires further research. 

9 Conclusions and Discussion 

In conclusion, we presented a brief review of recent developments for the FR/CPR schemes. Basic 
description and stability of the approach were discussed. Representative numerical examples were shown. 
Some key pacing items were mentioned. The FR/CPR framework appears to be promising and capable of 
contributing toward the goal of faster and higher fidelity CFD capabilities for more accurate flow field 
predictions. 
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